 |
 |
 |
|---|
Dr. Mehran Ghasabeh & Prof. Dr. Thomas Nagel
Chair of Soil Mechanics and Foundation Engineering
Institute of Geotechnics
TU Bergakademie Freiberg
Introduction¶
In nuclear waste disposal, a primary concern is the possibility of repository-induced or natural disturbances creating permeable flow paths within the initially impermeable host rock barrier. If such a flow path is to form and remain permanently open, it could serve as a transport pathway for radionuclides released from a waste package [1]. In geological nuclear waste disposal within argillaceous clay, fault activation can result from repository-scale thermally induced stress changes and a phenomenon known as thermal pressurization. Thermal pressurization refers to the increase in pore pressure caused by the thermal expansion of fluids trapped within the pores of low-permeability host rock. Without careful management through appropriate repository design, elevated repository temperatures and the resulting thermal pressurization could reach levels capable of inducing hydraulic fracturing or triggering shear activation of pre-existing fractures and minor faults [2, 3]. Pore pressure changes–caused by phenomena such as thermal pressurization, as demonstrated by Rutqvist et al. [4], or gas production [5, 6, 7]—can potentially reactivate faults, influencing mechanical responses and leading to associated changes in fault permeability.
This study presents the quantitative findings of a hydro-elasto-plastic model designed to examine the interactions between fault reactivation in clay host rock and the potential for enhanced fluid movement through a low-permeability medium. It addresses key challenges associated with the formation of permeable flow paths that could facilitate contaminant transport in otherwise low-permeability argillaceous rocks. These flow paths may develop due to fault reactivation of varying magnitudes, driven by thermal, hydraulic, and mechanical disturbances occurring during both operational and post-closure phases. To investigate these concerns, a permeability model incorporating embedded fractures [8] is employed. Rock masses inherently exhibit anisotropy in shear strength due to the presence of internal discontinuities. External and internal factors, such as pressure fluctuations, temperature changes, or mass removal, can significantly alter the effective stress state acting on these discontinuities. Consequently, plastic failure can occur even under compressive loads. Furthermore, shear-induced fracture dilation can lead to substantial changes in the hydraulic transmissivity tensor and fluid pressure within the rock mass. The primary objective of this study is to examine potential failure mechanisms in rock masses containing mechanical discontinuities such as faults, fractures, bedding planes, or other weak planar structures. Failure in such systems typically manifests in one of two ways: sliding failure along a weak plane or intrinsic failure of the rock mass itself. The rock matrix is assumed to exhibit either elastic behaviour or brittle failure, as described in [9].
We employ an elasto-plastic Coulomb failure-based weak plane constitutive model to systematically investigate shear failure mechanisms in rock masses containing mechanical discontinuities. This model explicitly accounts for the orientation and distribution of bedding features, such as embedded weak planes, to accurately simulate the rock's mechanical response. By incorporating these anisotropic characteristics, it enables a detailed analysis of stress redistribution, deformation, and failure within the rock mass.
The model is particularly suited for scenarios where a single fault serves as the primary mechanical discontinuity, offering valuable insights into fault behavior under varying loading conditions and its impact on rock stability and permeability. This approach provides a robust framework for understanding the influence of weak planes and fault orientation on shear failure, supporting the development of safer and more effective strategies for geological nuclear waste disposal."
In this study, a three-dimensional model domain with side lengths of $20$ meters is defined, incorporating a discontinuity dipping at an angle of $65$° (see Fig. 1). The fault is initially considered impermeable, consistent with field investigations conducted at the site. However, in the FM1 benchmark, the fault is assumed to be open and permeable within a $1$-meter radius around the well. Consequently, in the FM1 model, two distinct sets of material parameters are assigned to the fault, whereas in the FM2 model, a single material set is used.
The initial fluid pressure is set to $p_0=0.5$ MPa, based on site-specific measurements from the Mont Terri Laboratory. The stress field, though simplified from detailed site-specific estimations, assumes that the maximum compressive principal stress is vertical with a magnitude of $77$ MPa. The intermediate principal stress (also the maximum horizontal stress) is oriented perpendicular to the fault strike, with a magnitude of $66$ MPa. The minimum horizontal stress, parallel to the fault strike, has a magnitude of $3.3$ MPa [10]. However, since this stress is aligned with the fault, it does not influence fault activation. The model assumes uniform pressure and stress fields and neglects the influence of gravity. The 3D model domain is constrained by roller supports along each edge.
In this benchmark study, shear-induced hydraulic aperture changes and associated variations in fracture transmissivity are examined using two different fault models, FM1 and FM2. FM1 is developed based on insights from fault activation experiments conducted at the Tournemire site in Southern France [11]. Despite its relative simplicity, FM1 successfully replicates the measured responses observed at Tournemire. According to FM1, shear or tensile failure induces an aperture change of $\Delta{b}_{\text{ic}}=28$ $\mu\text{m}$, representing the creation of an open fault surface. This value is derived from modeling and analysis of similar field experiments at the Tournemire site. In FM1, elastic (normal) and plastic (shear) fault opening occurs only upon shear failure. In contrast, in FM2, an initial fracture aperture exists, allowing elastic (normal) opening from the outset with propagating pore pressure, while shear (dilatant) opening occurs only when shear failure is initiated.
The detailed model geometry including the locations of monitoring and anchors points is demonstrated in Fig. 2.
Description of Material Properties¶
This section discusses thoroughly how to calculate the mechanical and hydraulic parameters. The host rock is considered an elastic medium, whereas the fault is regarded as a weak plane embedded in elastoplastic media. In this study, the hydraulic conductivity of rock cracks is examined using the parallel plate model with cubic law. In the model, the fracture is assumed to be confined by two smooth, parallel walls separated by an aperture $b$. The flow is supposed to occur in a given fracture, with the flow rate directly proportional to the fluid pressure gradient and gravitational force, and transversely proportional to the dynamic viscosity $\mu$.
Mechanical Properties¶
The Young's modulus and the Poisson's ratio of the host rock matrix is calculated regarding the constitutive parameters, given in Table 1 by
$$ \tag{1} \begin{equation} E_{\text{M}} = \frac{9\, K\, G}{3\,K + G} \quad \text{and} \quad \nu_{\text{M}} = \frac{3\,K - 2\,\mu}{2(3\,K + \mu)}, \end{equation} $$
as $E_{\text{M}} = 6.1065$ GPa and $\nu_{\text{M}} = 0.3275$
Normal stiffness is evaluated considering the thickness of the fault and its Young's modulus accroding to the relations derived in [12]:
$$ \tag{2} \begin{equation} K_n = \frac{E\,(1- \nu)}{(1+\nu)(1-2\nu)}\frac{1}{t_{\text{el}}} \end{equation} $$
where the mean fracture distance is $t_{\text{el}}=0.021117$ m. Then, the Young's modulus of the fault is calculated as $E := 281.6$ kPa. Considering that the normal and shear stiffness of the fault are both equal, the Poisson's ratio is determined as $\nu = -0.5$.
Permeability Model¶
The fundamental concept of the model revealed here is the suitable depiction of individual fractures integrated into a continuous finite element. A finite element made up of a rock matrix and a number of $n$ fractures is seen in Figure 2. The breadth of each fracture, $a$, a characteristic size of the material, and the element size, $s$, determine the number of fractures in an element.
The fluid flow in the fracture is governed by Darcy's law, with a cubic dependency between the flow rate and the hydraulic conducting aperture. The fluid flow per fracture cross-section $q$ can be expressed as
$$ \tag{3} \begin{equation} q_i = -\frac{b^2}{12\,\mu} \left(p_{,i} - \rho_{\text{f}} g_i\right). \end{equation} $$
Consider first the flow phenomena through a single fracture. Liquid and gas flow will be calculated using Darcy’s law. The most important parameter in this law is the intrinsic permeability which can be calculated, assuming laminar flow, as
$$ \tag{4} \begin{equation} k_{\text{F}} = \frac{b^3}{12\,t_{\text{el}}}, \end{equation} $$
where $t_{\text{el}}$ indicates the thickness of the fault-embedded element that scales the permeability to correctly generate the cubic flow rule.
The transmisibility $T$ is the determined by
$$ \tag{5} \begin{equation} T_{\text{F}} = k\,s \quad \text{with} \quad k = k_{\text{M}} + \frac{b}{a} k_{\text{F}}. \end{equation} $$
The total permeability $k$ of the embedded fracture zone consists of the intrinsic matrix permeability $k_\text{M}$, the ratio of $b$ and the mean fracture distance $a$ as well as the fracture permeability $k_\text{F}$. $k_\text{M}$ is assumed to be nearly impermeable for the current benchmark example. In additon, there applies $s=a$ due to the existance of only one fracture in the benchmark. The expression of $T_{\text{F}} $ simplifies therefore to
$$ \tag{6} \begin{equation} T_{\text{F}} \approx \frac{b^2}{12}\,b. \end{equation} $$
When a set of $n$ fractures is included in a finite element, the equivalent intrinsic permeability of the element can be calculated as
$$ \tag{7} \begin{equation} k = k_{\text{M}} \boldsymbol{M} + \left( \frac{s-nb}{s}k_{\text{M}} + \frac{nb}{s}k_{\text{F}} \right) (\boldsymbol{I} - \boldsymbol{M} ) \quad \text{with} \quad n = \frac{s}{a}. \end{equation} $$
By replacing $k_{\text{F}}$ and $n$ in (7), we obtain the permeability as
$$ \tag{8} \begin{equation} k = k_{\text{M}} \boldsymbol{I} + \frac{b}{a} \left(\frac{b^2}{12} - k_{\text{M}} \right) (\boldsymbol{I} - \boldsymbol{M}), \end{equation} $$
where the hydraulic conducting aperture (opening) is additively decoupled into the initial aperture $b_{\text{i}}$, an induced fracture creation aperture change, $\Delta b_{\text{c}}$, an induced elastic aperture change, $\Delta b_{\text{e}}$, an induced shear dilation aperture change, $\Delta b_{\text{s}}$ as
$$ \tag{9} \begin{equation} b = b_{\text{i}} + \Delta b_{\text{e}} + \Delta b_{\text{ic}} . \end{equation} $$
In the current study, two different fault models, simply called FM1 and FM2, with different parameter sets are applied for dealing with shear-induced hydraulic aperture changes, and associated changes in fracture transmissivity. The models explain the different responses with respect to the evolution of the pressure along the fault path. The FM1 model is applied by introducing an initial creation aperture of $\Delta b_{\text{ic}} = 28 \mu$ m that comes from the modeling and evaluation of similar field experiments at the Tournemire site [5]. This parameter indicates a pre-set facture in the FM1 model, which is added as soon as the plastic deformation occurs. The induced elastic aperture change, $\Delta b_{\text{e}}$ is calculated by
$$ \tag{10} \begin{equation} \Delta b_{\text{e}} = a \, \langle \varepsilon_{\text{n}} - \varepsilon_{0} \rangle, \quad \text{with} \quad \varepsilon_{\text{n}}:=\boldsymbol{\varepsilon}:\boldsymbol{M} \end{equation} $$
where $\varepsilon_{0}$ indicates the threshold strain, and $\langle \varepsilon_{\text{n}} - \varepsilon_{0} \rangle$ states that the corresponding aperature change occurs only when $\varepsilon_{\text{n}} \gt \varepsilon_{0} $. The threshold strain $\varepsilon_{0}$ is determined through
$$ \tag{11} \begin{equation} b(t=0)=b_\text{i} + a \, \langle \varepsilon_\text{n}-\varepsilon_\text{0} \rangle, \end{equation} $$
while the normal strain $\varepsilon_\text{n}$ and the initial aperture $b_\text{i}$ are initially postulated to be zero. Therefore, the threshold strain can be determined by
$$ \tag{12} \begin{equation} \varepsilon_\text{0}= -\frac{b(t=0)}{a}. \end{equation} $$
In FM1, with $b(t=0) = 28 \mu\text{m}$, the threshold strain for the initially open fracture region is calculated as $\varepsilon_\text{0}:= -\cfrac{28^{-6}}{0.021117} = -1.32595\cdot 10^{-3}$, however for the fault is $\varepsilon_\text{0}:= 0$. In FM2, with $b(t=0) = 10 \mu\text{m}$, the threshold strain is calculated as $\varepsilon_\text{0}:= -\cfrac{10^{-5}}{0.021117} = -4.7355 \cdot 10^{-4}$.
Table 1
Constitutive parameters for modelling the induced slip of a fault
| Material | Parameter | Value (FM1) | Value (FM2) | Unit |
|---|---|---|---|---|
| Elastic parameters | Bulk modulus, $K$ | $5.9$ | $5.9$ | GPa |
| Shear modulus, $G$ | $2.3$ | $2.3$ | GPa | |
| Bulk density, $\rho_{\text{R}}$ | $2450$ | $2450$ | kg/m $^3$ | |
| Permeability, $k$ | 0 | 0 | m $^{2}$ | |
| Elastoplastic parameters | Normal Stiffness, $k_\text{n}$ | $20$ | $20$ | GPa/m |
| Shear stiffness, $k_\text{s}$ | $20$ | $20$ | GPa/m | |
| Weak plane cohesion, $c_{\text{wp}}$ | $0$ | $0$ | kPa | |
| Weak plane frictiona angle, $\phi_{\text{wp}}$ | $22$ | $22$ | $^\circ$ | |
| Weak plane dilatancy angle, $\psi_{\text{wp}}$ | $0$ | $10$ | $^\circ$ | |
| Initial creation aperture, $b_\text{i}$ | $28$ | $0$ | $\mu$ m | |
| Initial aperture, $b_0$ | $0$ | $10$ | $\mu$ m | |
| Hydraulic parameters | Permeability of intact rock, $k_\text{M}$ | $10^{-22}$ | $10^{-22}$ | $\text{m}^2$ |
| Initial aperature, $b_0$ | $10$ | $10$ | $\mu$ m | |
| Biot's coefficient, $\alpha$ | $1$ | $1$ | $-$ | |
| Fluid | Density, $\rho_{\text{f}}$ | $1000$ | $1000$ | kg/m $^3$ |
| Viscosity, $\mu$ | $1000$ | $1000$ | Pa $\cdot$ s |
Elasto-plastic weak plane models: Coulomb Failure Criterion¶
Rock masses can exhibit anisotropic shear strength due to the presence of one or more weakness planes. Variations in external or internal loads—such as changes in pressure, temperature, or mass removal—can alter the stress state acting on these planes. Consequently, plastic failure may occur even under compressive loading.
Failure can occur in two ways: (1) sliding along the weak plane or (2) intrinsic failure of the rock mass. The rock matrix is generally expected to behave elastically or, in extreme cases, in a brittle manner. In contrast, sliding failure is evaluated using the Coulomb failure criterion, applied to a specifically defined plane. This approach accounts for the presence of bedding features, commonly referred to as weak planes, within a three-dimensional framework. The behavior of jointed rock (including jointed rock masses, as considered in this paper) is anisotropic, non-linear, and dependent on the stress path due to the presence of discontinuities. Compared to intact rock, jointed rock generally exhibits lower stiffness, higher permeability, and reduced strength [12].
Along these planes, failure conditions are governed by the Coulomb failure envelope, which defines the transition from an intact state to failure based on the interaction between normal and shear stresses (see Fig. 3). Incorporating weak planes into the model enhances the accuracy of mechanical behavior simulations, particularly under stress conditions that may induce fault activation or shear deformation. This methodology is crucial for geomechanical analyses, including fault reactivation, fracture propagation, and subsurface engineering applications.
The stress point's position on the failure envelope is determined using a non-associated flow rule for shear failure. The weak-plane orientation is specified by the Cartesian components of a unit normal vector in the global $x$-, $y$-, and $z$- axes. A local reference coordinate system is defined with basis vectors $e0§, $e1$, and e3e3. In this formulation, the normal vector $e_0$ and the tangential vector e1e1 are computed with respect to a fault-embedded medium dipping at $65^\circ$. Thus, in the local coordinate system, the basis vectors are $e_0=(−0.906,0.423)$ and $e_1=(−0.423,−0.906)$.
The Coulomb yield and potential functions,
$$ \tag{15} \begin{equation} F^{\mathrm{wp}} = \sigma^{\mathrm{wp}}_n \tan \phi^{\mathrm{wp}} - c^{\mathrm{wp}} + \tau^{\mathrm{wp}} \end{equation} $$
and
$$ \tag{16} \begin{equation} G^{\mathrm{wp}}_F = \sigma^{\mathrm{wp}}_n \tan \phi^{\mathrm{wp}} - c^{\mathrm{wp}} + \tau^{\mathrm{wp}} \end{equation} $$
are used to extend a multi-surface plasticity envolving weakness plane with normal $n^{\mathrm{wp}}$, regarding non-associated Coulomb behaviour. In above equations, the normal and the shear stresses on the fault path are calculated by
$$ \begin{equation} \tag{17} \sigma^{\mathrm{wp}}_{\mathrm{n}} = \boldsymbol{\sigma} : \boldsymbol{n} \otimes \boldsymbol{n} \end{equation} $$
and
$$ \tag{18} \begin{equation} \tau^{\mathrm{wp}}= \sqrt{\sigma^{\mathrm{T}}\sigma:\boldsymbol{n} \otimes \boldsymbol{n}- (\boldsymbol{\sigma} : \boldsymbol{n} \otimes \boldsymbol{n})^2} \end{equation} $$
Implementation of Constitutive Relations¶
The open-source software packages OpenGeoSys [15], MFront [16] and MFrontGenericInterfaceSupport (MGIS) [17] are among the systematic techniques we identify as being able to set up such modeling workflows. In this work, the elasto-plasticity approach with an embedded weakness plane for hydro-mechanical fault reactivation modeling is implemented in MFront. The material knowledge platform MFront and the generic interface library MGIS were developed to make the process of constitutive model implementation less error-prone as well as more reproducible, efficient, and user-friendly across a wide range of process codes. MFront features a code generator that converts a collection of closely related domain-specific languages into plain C++ sources. It makes coding and integration of constitutive equations easier to implement as well as portable between different solvers, thus contributing to the creation of reproducible simulations and efficient code. With the help of MGIS, MFront can be integrated in a large range of solvers. The OpenGeoSys project, coupled with MFront, can provide an efficient and robust numerical tool for geoscientific applications that simulate many multiphysics issues.
Calculation of the normal and the shear displacements¶
The normal and the shear displacements on the lines located on the right-hand, middle, and left-hand side of the fracture-embedded fault, see Fig. 3, are calculated by
$$ \begin{equation} u_n = {u}_x \cos \theta + {u}_y \sin \theta \quad \text{and} \quad u_s =-{u}_x \sin \theta + {u}_y \cos \theta \end{equation} $$
As represented in Fig. 3, the normal and the shear displacements on the Line AB, Line CD, Line EF, and the points M and N are calculated withe respect to the normal vector, $e_0 := (\cos(155), \sin(155), 0.) = (-0.9063, 0.4226, 0.)$
Installation of OpenGeoSys with PETSc for Parallel Computing¶
Then get the source code by cloning the source code repository, and the specific feature branch with Git, using the following command:
git clone -b FaultSlipTestFM1 --single-branch https://gitlab.opengeosys.org/MehranGhasabeh/ogs-mg.git
For more information, please refer to the official documentation at https://www.opengeosys.org/docs/devguide/getting-started/get-the-source-code/.
To install OpenGeoSys with PETSc is then required to configure and finally build the configuration. For configuration, run cmake with preset in the source code directory:
cmake --preset release-petsc -DBUILD_SHARED_LIBS=ON -DOGS_BUILD_HDF5=ON -DOGS_USE_MKL=ON -DOGS_EIGEN_PARALLEL_BACKEND=OpenMP -DOGS_USE_MFRONT=ON -DCPM_SOURCE_CACHE=$HOME/.cache/CPM -DOGS_PETSC_CONFIG_OPTIONS="--with-debugging=0 --download-superlu --download-superlu_dist --download-scalapack --download-mumps
For more information, please refer to the official documentation at
To build with the ninja tool on the shell go to your previously configured build directory and type ninja.
For detailed information please refer to the the official documentation at https://www.opengeosys.org/docs/devguide/getting-started/build/.
Mesh Generation¶
geo3D.geo defines the geometry, assigns physical groups,
and applies a refinement field to create a bias in the mesh density around the fault center.
import numpy as np
# Parameters for the line
original_length = -15.0 # original total length of the line in meters
new_length = -10.0 # desired total length after scaling
num_elements = 18 # number of elements
bias_factor = 1.25 # growth factor for bias (applied to the scaled line only)
# Generate original node coordinates for a uniform 3-meter line (no bias)
original_node_coords = np.linspace(0, original_length, num_elements + 1)
# Calculate the length of the first element for the biased 2-meter line
l0 = new_length * (bias_factor - 1) / (bias_factor**num_elements - 1)
# Generate biased node coordinates for the new 2-meter line
scaled_node_coords = [0.0]
for i in range(num_elements):
scaled_node_coords.append(scaled_node_coords[-1] + l0 * bias_factor**i)
def format_to_16_chars(number):
# Calculate the integer part length
integer_part_length = len(str(int(abs(number))))
# Calculate the total width by subtracting the integer part length (if > 0)
total_width = 18 - integer_part_length if integer_part_length > 0 else 18
# Calculate the number of decimal places to keep the formatted number within the total width
decimal_places = (
total_width - 1
) # Subtract 1 for the decimal point and 1 extra for formatting
decimal_places = max(0, decimal_places) # Ensure non-negative decimal places
# Format the number with the calculated decimal places
formatted_value = f"{number:>{18}.{decimal_places}f}"
return formatted_value
# Print each pair in the format "old coordinate : new coordinate"
for old, new in zip(original_node_coords, scaled_node_coords):
print(f"{format_to_16_chars(old)} : {format_to_16_chars(new)},")
1_extrusion.py script is utilized to refine the mesh specifically around the injection point located on the fault. It leverages the biased meshing generated by the previous code block, ensuring that the mesh density is higher in critical regions near the injection point.
msh file generated by Gmsh to OGS mesh format, generate quadratic mesh and create the required geometrical mesh information including the definition of regions with specific material IDs and the identification of boundaries.
import subprocess
import os
# Create a directory for saving VTU files
output_dir = "inputMeshes"
os.makedirs(output_dir, exist_ok=True)
# Remove unnecessary files
commands = [
"rm -f \\#* *~ *.msh",
"gmsh -3 ./GEOs/geo3D.geo -o ./GEOs/geo3D.msh -format msh2",
"python3 1_extrusion.py",
"~/build/release-petsc/bin/GMSH2OGS -i ./GEOs/geo3D.msh -o ./GEOs/geo.vtu -e --gmsh2_physical_id",
"~/build/release-petsc/bin/NodeReordering -i ./GEOs/geo.vtu -o ./GEOs/geo_reorder.vtu",
"~/build/release-petsc/bin/editMaterialID -i ./GEOs/geo_reorder.vtu -o ./GEOs/geo_reorder.vtu -r -m 0 -n 0",
"~/build/release-petsc/bin/editMaterialID -i ./GEOs/geo_reorder.vtu -o ./GEOs/geo_reorder.vtu -r -m 1 -n 0",
"~/build/release-petsc/bin/editMaterialID -i ./GEOs/geo_reorder.vtu -o ./GEOs/geo_reorder.vtu -r -m 2 -n 1",
"~/build/release-petsc/bin/editMaterialID -i ./GEOs/geo_reorder.vtu -o ./GEOs/geo_reorder.vtu -r -m 3 -n 0",
"~/build/release-petsc/bin/editMaterialID -i ./GEOs/geo_reorder.vtu -o ./GEOs/geo_reorder.vtu -r -m 4 -n 0",
f"mv ./GEOs/geo_reorder.vtu {output_dir}/geo_domain_3D.vtu",
f"~/build/release-petsc/bin/createQuadraticMesh -i {output_dir}/geo_domain_3D.vtu -o {output_dir}/geo_domain_3D_q8.vtu",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_top.vtu -x 0. -y -1. -z 0.",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_bottom.vtu -x 0. -y 1. -z 0.",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_back.vtu -x 1. -y 0. -z 0.",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_front.vtu -x -1. -y 0. -z 0.",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_left.vtu -x 0. -y 0. -z -1.",
f"~/build/release-petsc/bin/ExtractSurface -i {output_dir}/geo_domain_3D_q8.vtu -o {output_dir}/geometry_right.vtu -x 0. -y 0. -z 1.",
"python3 materialIDs3D.py",
]
# Execute commands
for cmd in commands:
print(f"Executing: {cmd}")
result = subprocess.run(cmd, shell=True, capture_output=True, text=True)
if result.returncode != 0:
print(f"Error executing command: {cmd}")
print(result.stderr)
else:
print(result.stdout)
import shutil
# Define the output directory
output_dir_FM1 = "./FM1/inputMeshes"
output_dir_FM2 = "./FM2/inputMeshes"
# Define the commands
cmd1_FM1 = f"~/build/release-petsc/bin/partmesh -i {output_dir_FM1}/geo_domain_3D_q8_FM1.vtu --ogs2metis"
cmd2_FM1 = f"~/build/release-petsc/bin/partmesh -n 24 -m -i {output_dir_FM1}/geo_domain_3D_q8_FM1.vtu -- {output_dir_FM1}/geometry_*.vtu"
cmd1_FM2 = f"~/build/release-petsc/bin/partmesh -i {output_dir_FM2}/geo_domain_3D_q8_FM2.vtu --ogs2metis"
cmd2_FM2 = f"~/build/release-petsc/bin/partmesh -n 24 -m -i {output_dir_FM2}/geo_domain_3D_q8_FM2.vtu -- {output_dir_FM2}/geometry_*.vtu"
# Run the commands
try:
print("Running first command...")
subprocess.run(cmd1_FM1, shell=True, check=True) # Run the first command
print("First command completed (FM1).")
print("Running second command...")
subprocess.run(cmd2_FM1, shell=True, check=True) # Run the second command
print("Second command completed (FM1).")
# Move generated files to output_dir
print("Organizing generated files...")
for file in os.listdir("."): # Current directory
if file.endswith((".bin", ".mesh")): # Adjust file types as needed
shutil.move(file, output_dir_FM1)
print("All generated files moved to the FM1 output directory.")
except subprocess.CalledProcessError as e:
print(f"An error occurred: {e}")
# Run the commands
try:
print("Running first command...")
subprocess.run(cmd1_FM2, shell=True, check=True) # Run the first command
print("First command completed (FM1).")
print("Running second command...")
subprocess.run(cmd2_FM2, shell=True, check=True) # Run the second command
print("Second command completed (FM1).")
# Move generated files to output_dir
print("Organizing generated files...")
for file in os.listdir("."): # Current directory
if file.endswith((".bin", ".mesh")): # Adjust file types as needed
shutil.move(file, output_dir_FM2)
print("All generated files moved to the FM1 output directory.")
except subprocess.CalledProcessError as e:
print(f"An error occurred: {e}")
To define the pressure injection boundary mesh, we use ParaView (https://www.paraview.org/). The process involves selecting and extracting the region where the pressure will be applied, followed by verifying its consistency with the bulk mesh using the identifySubdomains tool from the OpenGeoSys framework. The detailed steps are as follows:
- Open the mesh file in ParaView.
- Navigate to
Filters > Alphabetical > Extract Selectionto select the region of interest for the pressure boundary. - After extracting the selection, save the resulting data using
ParaView > File > Save Data. - Save the file in
.vtuformat to ensure compatibility with the next steps. - Utilize the
identifySubdomainstool to ensure that the extracted subdomain meshes are part of the bulk mesh and to map their corresponding node and element IDs:~/build/release-petsc/bin/identifySubdomains -m geo_domain_3D_q8_FM1.vtu geometry_pinj.vtuThis command checks the consistency between the bulk mesh (geo_domain_3D_q8_FM1.vtu) and the subdomain (geometry_pinj.vtu), and generates files that contain:
bulk_node_ids: Nodes in the subdomain mapped to the bulk mesh.bulk_element_ids: Elements in the subdomain mapped to the bulk mesh. For additional details on using the identifySubdomains tool, refer to the official OpenGeoSys documentation at https://www.opengeosys.org/docs/tools/meshing-submeshes/identifysubdomains/.
Project File Structure Overview¶
This section presents a structured overview of the project file hierarchy, outlining the organization of directories and files to enhance clarity, streamline navigation, and facilitate a comprehensive understanding of the model implementation.
Mesh inputs¶
Processes definition in an OGS .prj File¶
This section provides a demonstration of how to define processes, including the constitutive relation, material properties, process variables, secondary variables, and initial stress within an OGS .prj file.
Material properties definition as parameters in an OGS .prj File¶
Definition of displacement BCs in an OGS .prj file¶
Definition of Pressure BCs in an OGS .prj file¶
Definition of permeability model inputs in an OGS .prj file¶
EmbeddedFracturePermeability model inputs around injection well in an OGS .prj file using FM1 formulation.
EmbeddedFracturePermeability model inputs for fault in an OGS .prj file using FM1 formulation.
EmbeddedFracturePermeability model inputs for fault in an OGS .prj file using FM2 formulation.
import vtuIO
import math
import csv
from pathlib import Path
import pandas as pd
import numpy as np
import pyvista as pv
import matplotlib.pyplot as plt
## To run the simulations through the jupyternotebook, the following lib is imported
from ogs6py.ogs import OGS
# Some plot settings
plt.rcParams["lines.linewidth"] = 3.0
plt.rcParams["lines.color"] = "black"
plt.rcParams["legend.frameon"] = True
plt.rcParams["font.family"] = "serif"
plt.rcParams["legend.fontsize"] = 18
plt.rcParams["font.size"] = 18
plt.rcParams["axes.spines.right"] = True
plt.rcParams["axes.spines.top"] = True
plt.rcParams["axes.spines.left"] = True
plt.rcParams["axes.spines.bottom"] = True
plt.rcParams["axes.axisbelow"] = True
FM1 Model¶
Injection pressure¶
dt = np.array(
[0, 1, 1, 22.3, 1, 76, 1, 56, 1, 55, 1, 53, 1, 49, 1, 102, 1, 32, 1, 47, 100, 200]
) # seconds
p = np.array(
[
0.5,
0.5,
0.7446,
0.7446,
1.919,
1.919,
3.627,
3.627,
4.094,
4.094,
4.511,
4.5411,
4.99,
4.99,
5.484,
5.484,
6.302,
6.302,
3.382,
3.382,
3.382,
3.382,
]
) # MPa
fig, ax = plt.subplots(figsize=(9, 5))
ax.plot(
dt.cumsum()[1:] - dt[1],
p[1:],
label="Injected pressure",
color="red",
)
ax.set_xlabel("$t-t_\\mathrm{ic}$ / s")
ax.set_ylabel("Well-Pressure / MPa")
ax.legend()
ax.grid(True)
fig.tight_layout()
Run simulation¶
To run the simulation through the jupyternotebook, at first it is required to call from ogs6py.ogs import OGS, then a directory for the output is created.
out_dir_FM1 = Path(os.environ.get("OGS_TESTRUNNER_OUT_DIR", "_out_FM1"))
if not out_dir_FM1.exists():
out_dir_FM1.mkdir(parents=True)
model_fm1=OGS(INPUT_FILE="FST_EPKd_FM1.prj", PROJECT_FILE="FST_EPKd_FM1.prj")
Following that, the input and project files are defined.
PATH_OGS = "~/build/release-petsc/bin"
model_fm1.run_model(path=PATH_OGS, logfile=str(out_dir_FM1 / "out.txt"), args=f"-o {out_dir_FM1}")
Import ENSI data from the Excel file¶
injectionPressure = pd.read_excel(
"Injectionpressure.xlsx", usecols="A,E", skiprows=1, header=None, engine="openpyxl"
)
# print(injectionPressure[4])
FM1 = pd.read_excel(
"D19_TaskB_Step1_output_ENSI_FM1.xlsx",
sheet_name="Time dependent output",
usecols="B:W",
skiprows=range(0, 12),
header=None,
engine="openpyxl",
)
FM2 = pd.read_excel(
"D19_TaskB_Step1_output_ENSI_FM2.xlsx",
sheet_name="Time dependent output",
usecols="B:W",
skiprows=range(0, 12),
header=None,
engine="openpyxl",
)
profiles_along_fault_dip = pd.read_excel(
"D19_TaskB_Step1_output_ENSI_FM1.xlsx",
sheet_name="Profiles along fault dip",
usecols="A:J",
skiprows=range(0, 3),
header=None,
engine="openpyxl",
)
profiles_along_fault_dip_FM2 = pd.read_excel(
"D19_TaskB_Step1_output_ENSI_FM2.xlsx",
sheet_name="Profiles along fault dip",
usecols="A:J",
skiprows=range(0, 3),
header=None,
engine="openpyxl",
)
Import the data from PVD files¶
pvdfile_elastop_Kd_FM1 = vtuIO.PVDIO("./_out_FM1/FST_EPKd_FM1.pvd", dim=3)
Coordinates of injection, anchors, and monitoring Points¶
Anchors at poinsts at $(0, 0, 0.25)$ and $(0, 0, -0.25)$, injection at poinst at $(0, 0, 0)$, and moniroting point at $(0, -0.63, -1.36)$
anchors_point_1 = (10.0, 10.25, 0.0)
anchors_point_2 = (10.0, 9.75, 0.0)
hostRock_ptu = (10.0, 15.0, 0.0)
hostRock_ptb = (10.0, 5.0, 0.0)
##
fault_thickness = 0.0233
deltaL = 1.5
dip_angle = 65.0
dx = (
fault_thickness
/ 2.0
* math.sin(math.radians(dip_angle))
* math.sin(math.radians(dip_angle))
)
dy = (
fault_thickness
/ 2.0
* math.sin(math.radians(dip_angle))
* math.cos(math.radians(dip_angle))
)
##
pinjection_point = (10.0, 10.0, 0.0)
pinjection_point_left = (pinjection_point[0] - dx, pinjection_point[1] + dy, 0.0)
pinjection_point_right = (pinjection_point[0] + dx, pinjection_point[1] - dy, 0.0)
monitoring_point_p3 = (
pinjection_point[0] - deltaL * math.cos(math.radians(dip_angle)),
pinjection_point[1] - deltaL * math.sin(math.radians(dip_angle)),
0.0,
)
monitoring_point_p3_left = (
monitoring_point_p3[0] - dx,
monitoring_point_p3[1] + dy,
0.0,
)
monitoring_point_p3_right = (
monitoring_point_p3[0] + dx,
monitoring_point_p3[1] - dy,
0.0,
)
monitoring_point_p2 = (10.0, 10.0, -1.5)
monitoring_point_p2_left = (
monitoring_point_p2[0] - dx,
monitoring_point_p2[1] + dy,
-1.5,
)
monitoring_point_p2_right = (
monitoring_point_p2[0] + dx,<code>
monitoring_point_p2[1] - dy,
-1.5,
)
flow_rate_pt_0 = (9.977722711419313, 9.977209706092465, 0.0)
flow_rate_pt_1 = (9.990305916447044, 10.00419447635663, 0.0)
flow_rate_pt_2 = (10.00288912147477, 10.03117924662077, 0.0)
flow_rate_pt_3 = (9.997110878523848, 9.968820753376264, 0.0)
flow_rate_pt_4 = (10.00969408355157, 9.995805523640405, 0.0)
flow_rate_pt_5 = (10.02227728857931, 10.02279029390458, 0.0)
flow_rate_pt_6 = (9.977722711419313, 9.977209706092465, -0.0458621757848043)
flow_rate_pt_7 = (9.990305916447044, 10.00419447635663, -0.0458621757848043)
flow_rate_pt_8 = (10.00288912147477, 10.03117924662077, -0.0458621757848043)
flow_rate_pt_9 = (9.997110878523848, 9.968820753376264, -0.0458621757848043)
flow_rate_pt_10 = (10.00969408355157, 9.995805523640405, -0.0458621757848043)
flow_rate_pt_11 = (10.02227728857931, 10.02279029390458, -0.0458621757848043)
flowRate_points = {
"pt0": flow_rate_pt_0,
"pt1": flow_rate_pt_1,
"pt2": flow_rate_pt_2,
"pt3": flow_rate_pt_3,
"pt4": flow_rate_pt_4,
"pt5": flow_rate_pt_5,
"pt6": flow_rate_pt_6,
"pt7": flow_rate_pt_7,
"pt8": flow_rate_pt_8,
"pt9": flow_rate_pt_9,
"pt10": flow_rate_pt_10,
"pt11": flow_rate_pt_11,
}
pressure_points = {
"pt0": pinjection_point,
"pt1": monitoring_point_p3,
"pt2": monitoring_point_p2,
}
displacement_points = {
"pt0": anchors_point_1,
"pt1": anchors_point_2,
"pt2": pinjection_point_left,
"pt3": pinjection_point_right,
"pt4": monitoring_point_p3_left,
"pt5": monitoring_point_p3_right,
"pt6": monitoring_point_p2_left,
"pt7": monitoring_point_p2_right,
"pt8": hostRock_ptu,
"pt9": hostRock_ptb,
}
stress_points = {
"pt0": pinjection_point,
"pt1": monitoring_point_p2,
"pt2": monitoring_point_p3,
"pt3": hostRock_ptu,
}
Calling the results array of the interpolated pressure at the monitoring points.¶
timeSteps = pvdfile_elastop_Kd_FM1.timesteps
results = {}
results["pressure"] = pvdfile_elastop_Kd_FM1.read_time_series(
"pressure_interpolated", pts=pressure_points
)
Save the pressure data in a csv file¶
# Output CSV file
output_csv = "pressure_results_FM1.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 3
# Write header row: "Time" + Node IDs (e.g., pt0, pt1, pt2)
header = ["Time"] + [
f"pt{node_id}" for node_id in range(node_ids)
] # Adjust as needed
writer.writerow(header)
# Write each row: time step + scaled pressures for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Scale the pressure value by 1e-6 and append it
scaled_pressure = results["pressure"][f"pt{node_id}"][i] * 1e-6
row.append(scaled_pressure)
writer.writerow(row)
Calling the results array of the mass flow rate at the injection points.¶
results["MassFlowRate"] = pvdfile_elastop_Kd_FM1.read_time_series(
"MassFlowRate", pts=flowRate_points
)
Save the mass flow rate data in a csv file¶
# Output CSV file
output_csv = "MassFlowRate_results_FM1.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 12
# Write header row: "Time" + Node IDs (e.g., pt0 ~ pt11)
header = ["Time"] + [
f"pt{node_id}" for node_id in range(node_ids)
] # Adjust as needed
writer.writerow(header)
# Write each row: time step + scaled pressures for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Scale the pressure value by 1e-6 and append it
scaled_pressure = results["MassFlowRate"][f"pt{node_id}"][i]
row.append(scaled_pressure)
writer.writerow(row)
Calling the results array of the displacements at the anchors.¶
results["displacement"] = pvdfile_elastop_Kd_FM1.read_time_series(
"displacement", pts=displacement_points
)
Save the displacement data in a csv file¶
# Output CSV file
output_csv = "displacement_results_FM1.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 10 # Number of nodes (pt0 ~ pt9)
# Write header row: "Time" + Displacement Components (e.g., pt0_ux, pt0_uy, pt0_uz)
header = ["Time"]
for node_id in range(node_ids):
header.extend([f"pt{node_id}_ux", f"pt{node_id}_uy", f"pt{node_id}_uz"])
writer.writerow(header)
# Write each row: time step + displacement components for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Extract the displacement vector for the node and split into components
displacement_vector = results["displacement"][f"pt{node_id}"][i]
ux, uy, uz = (
displacement_vector[0],
displacement_vector[1],
displacement_vector[2],
)
row.extend([ux, uy, uz])
writer.writerow(row)
results["stress"] = pvdfile_elastop_Kd_FM1.read_time_series("sigma", pts=stress_points)
Save the stress data in a csv file¶
# Output CSV file
output_csv = "stress_results_FM1.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 5 # Number of nodes (pt0 ~ pt9)
# Write header row: "Time" + Displacement Components (e.g., pt0_ux, pt0_uy, pt0_uz)
header = ["Time"]
for node_id in range(node_ids):
header.extend([f"pt{node_id}_sigxx", f"pt{node_id}_sigyy", f"pt{node_id}_sigzz", f"pt{node_id}_sigxy", f"pt{node_id}_sigyz", f"pt{node_id}_sigxz"])
writer.writerow(header)
# Write each row: time step + displacement components for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Extract the displacement vector for the node and split into components
stress_tensor = results["stress"][f"pt{node_id}"][i]
sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_xz = (
stress_tensor[0],
stress_tensor[1],
stress_tensor[2],
stress_tensor[3],
stress_tensor[4],
stress_tensor[5]
)
row.extend([sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_xz])
writer.writerow(row)
Modeling results for FM1 benchmark¶
Pressure at P3: $(0, -0.63, -1.36)$¶
fig, ax = plt.subplots(figsize=(14, 6.5))
# Read the CSV file
csv_file = "pressure_results_FM1.csv"
pressureData = pd.read_csv(csv_file)
ax.plot(
FM1[1],
FM1[4] * 1e-6,
label="pressure at monitoring point, ENSI",
color="blue",
ls="--",
marker="o",
)
ax.plot(
pressureData["Time"],
pressureData["pt1"],
# instead the following two lines are used to plot the results
# timeSteps,
# results["pressure"]["pt1"] * 1e-6,
label="Pressure at monitoring Point, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Pressure at monitoring point P3 ($0.$, $-0.63$, $-1.36$)")
ax.set_xlabel("Time / s")
ax.set_ylabel("Pressure / MPa")
ax.legend()
ax.grid(True)
fig.tight_layout()
fig, ax = plt.subplots(figsize=(14, 6.5))
# Read the CSV file
csv_file = "pressure_results_FM1.csv"
pressureData = pd.read_csv(csv_file)
ax.plot(
FM1[1],
FM1[3] * 1e-6,
label="pressure at monitoring point, ENSI",
color="blue",
ls="--",
marker="o",
)
ax.plot(
pressureData["Time"],
pressureData["pt2"],
# instead the following two lines are used to plot the results
# timeSteps,
# results["pressure"]["pt1"] * 1e-6,
label="Pressure at monitoring Point, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Pressure at monitoring point P2 ($-1.5$, $0.$, $0.$)")
ax.set_xlabel("Time / s")
ax.set_ylabel("Pressure / MPa")
ax.legend()
ax.grid(True)
fig.tight_layout()
Pressure along fault dip at $t= 317$¶
start_point_left = (5.32527, 0.0, 0.0)
end_point_left = (14.6514, 20.0, 0.0)
start_point_right = (5.34857, 0.0, 0.0)
end_point_right = (14.6747, 20.0, 0.0)
start_point_mid = (
start_point_left[0] + (start_point_right[0] - start_point_left[0]) / 2.0,
0.0, 0.0
)
end_point_mid = (
end_point_left[0] + (end_point_right[0] - end_point_left[0]) / 2.0,
20.0, 0.0
)
stop = np.sqrt(
(end_point_mid[0] - start_point_mid[0]) ** 2.0
+ (end_point_mid[1] - start_point_mid[1]) ** 2.0
)
def divide_line(start_point, end_point, n):
fault_axis = []
for i in range(0, n + 1):
x = start_point[0] + (end_point[0] - start_point[0]) * i / n
y = start_point[1] + (end_point[1] - start_point[1]) * i / n
fault_axis.append([x, y, 0.0])
return fault_axis
npts = 1000
fpoints_mid = divide_line(start_point_mid, end_point_mid, npts)
fpoints_left = divide_line(start_point_left, end_point_left, npts)
fpoints_right = divide_line(start_point_right, end_point_right, npts)
rad = np.linspace(start=0, stop=stop, num=npts + 1)
import matplotlib.pyplot as plt
# Create the plot figure and axis
fig, ax = plt.subplots(figsize=(14, 6.5))
# Plot the first dataset (ENSI pressure along the fault dip) after filtering out zeros
nonzero_indices_1 = profiles_along_fault_dip[1] != 0 # Find non-zero indices
ax.plot(
profiles_along_fault_dip[0][nonzero_indices_1], # Filter x-values
profiles_along_fault_dip[1][nonzero_indices_1]
* 1e-6, # Filter y-values and convert to MPa
label="Pressure along the fault dip at $t=317$ s, ENSI",
color="blue",
ls="--",
marker="o",
)
# Define the time for the second dataset
itime = 317
# Plot the second dataset (TUBAF pressure along the fault tip) after filtering out zeros
rad_data = rad[::5]
pressure_data = (
pvdfile_elastop_Kd_FM1.read_set_data(
itime, "pressure_interpolated", pointsetarray=fpoints_mid
)[::5]
* 1e-6
)
nonzero_indices_2 = pressure_data != 0 # Find non-zero indices
ax.plot(
rad_data[nonzero_indices_2], # Filter x-values
pressure_data[nonzero_indices_2], # Filter y-values
label=f"Pressure along the fault tip at $t = {itime}$ s, TUBAF",
color="red",
ls="-",
marker="d",
)
# Set plot title and axis labels
ax.set_title("Pressure along the fault dip at $t=317$ s")
ax.set_xlabel("Fault dip ($r$) / m")
ax.set_ylabel("Pressure / MPa")
# Add legend and grid
ax.legend()
ax.grid(True)
# Adjust layout
fig.tight_layout()
# Show the plot
plt.show()
Injection Rate¶
fig, ax = plt.subplots(figsize=(14.5, 6.5))
# Read the CSV file
csv_file = "MassFlowRate_results_FM1.csv"
MassFlowRateData = pd.read_csv(csv_file)
ax.plot(
FM1[1],
FM1[2],
label="Injection rate, ENSI",
color="blue",
ls="--",
marker="o",
)
# Initialize the total mass flow rate
total_mass_flow_rate = 0
# Loop through keys pt0 to pt11 and sum the values
for i in range(12): # From pt0 to pt11
key = f"pt{i}"
total_mass_flow_rate += results["MassFlowRate"][key] * 60.0 * 2.0
MassFlowRateData["TotalMassFlowRate"] = (
MassFlowRateData.iloc[:, 1:].sum(axis=1) * 60.0 * 2.0
) # Sum across all columns except "Time"
ax.plot(
MassFlowRateData["Time"],
MassFlowRateData["TotalMassFlowRate"],
# timeSteps,
# total_mass_flow_rate,
label="Injection rate, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Injection rate versus Time")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Injection Rate [L / min]")
ax.legend(loc="lower right")
ax.grid(True)
fig.tight_layout()
Differential displacements between anchors¶
fig, ax1 = plt.subplots(figsize=(14.5, 6.5))
csv_file = "displacement_results_FM1.csv"
displacementData = pd.read_csv(csv_file)
ax1.plot(
FM1[1],
FM1[6] * 1000.0,
label="Differential displacement ($dy$), ENSI",
color="blue",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM1[1],
FM1[7] * 1000.0,
label="Differential displacement ($dz$), ENSI",
color="red",
ls="--",
marker="o",
)
# Line plot with markers
ax1.plot(
displacementData["Time"],
(displacementData["pt0_ux"] - displacementData["pt1_ux"]) * 1e3,
# timeSteps,
# (results["displacement"]["pt0"].T[0]
# - results["displacement"]["pt1"].T[0]) * 1e3,
label="Differential displacement ($dy$), TUBAF",
color="magenta",
ls="-",
marker="D",
)
ax1.plot(
displacementData["Time"],
(displacementData["pt0_uy"] - displacementData["pt1_uy"]) * 1e3,
# timeSteps,
# (results["displacement"]["pt0"].T[1]
# - results["displacement"]["pt1"].T[1]) * 1e3,
label="Differential displacement ($dz$), TUBAF",
color="green",
ls="-",
marker="D",
)
ax1.set_title("Relative displacement of anchors")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Rotation matrix with respect to the fault dip angle¶
theta = 155
rotation_matrix = np.array(
[
[math.cos(math.radians(theta)), math.sin(math.radians(theta))],
[-math.sin(math.radians(theta)), math.cos(math.radians(theta))],
]
)
print(rotation_matrix)
[[-0.90630779 0.42261826] [-0.42261826 -0.90630779]]
Calculation of the normal and shear displacment along the fault dip¶
def displacement_ns(file, theta, pt_left, pt_right):
u_nd_left = (
file["displacement"][pt_left].T[0] * math.cos(math.radians(theta))
+ file["displacement"][pt_left].T[1] * math.sin(math.radians(theta))
)
u_sd_left = -file["displacement"][pt_left].T[0] * math.sin(
math.radians(theta)
) + file["displacement"][pt_left].T[1] * math.cos(math.radians(theta))
##
u_nd_right = file["displacement"][pt_right].T[0] * math.cos(
math.radians(theta)
) + file["displacement"][pt_right].T[1] * math.sin(math.radians(theta))
u_sd_right = -file["displacement"][pt_right].T[0] * math.sin(
math.radians(theta)
) + file["displacement"][pt_right].T[1] * math.cos(math.radians(theta))
return abs(u_nd_left) + abs(u_nd_right), abs(u_sd_left) + abs(u_sd_right)
u_nd, u_sd = displacement_ns(results, 155, "pt2", "pt3")
Fault displacement at P1: $(0, 0, 0)$¶
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[8] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM1[1],
FM1[10] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="-",
marker="d",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="-",
marker="d",
)
ax1.set_title("Fault displacement at $(0, 0, 0)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
u_nd, u_sd = displacement_ns(results, 155, "pt4", "pt5")
Fault displacement at P3: $(0., -0.63, -1.36)$¶
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[14] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM1[1],
FM1[16] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_title("Fault displacement at P2 $(0., -0.63, -1.36)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
u_nd, u_sd = displacement_ns(results, 155, "pt6", "pt7")
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[11] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM1[1],
FM1[13] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_title("Fault displacement at P2 $(-1.5, 0., 0.)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Stress in host rock at $(15. , 10., 0.)$¶
fig, ax1 = plt.subplots(2, 2, figsize=(20, 9))
csv_file = "stress_results_FM1.csv"
stressData = pd.read_csv(csv_file)
# Top-left plot
ax1[0, 0].plot(
stressData["Time"],
stressData["pt3_sigxx"],
#timeSteps,
#results["stress"]["pt3"].T[0] * 1e-6,
label="$\\sigma_{xx}$, TUBAF",
color="green",
ls="--",
marker="o",
)
ax1[0, 0].set_title("$\\sigma_{xx}$ vs Time at $(0., 5., 0.)$")
ax1[0, 0].set_xlabel("Time / s")
ax1[0, 0].set_ylabel("$\\sigma_{xx}$ / MPa")
ax1[0, 0].legend()
ax1[0, 0].grid(True)
# Create a second y-axis for the second data series
ax2_00 = ax1[0, 0].twinx()
color = "tab:red"
ax2_00.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2_00.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2_00.tick_params(axis="y", labelcolor=color)
ax2_00.set_xlim([0, 850])
# Top-right plot
ax1[0, 1].plot(
stressData["Time"],
stressData["pt3_sigxy"],
#timeSteps,
#results["stress"]["pt3"].T[1] * 1e-6,
label="$\\sigma_{yy}$, TUBAF",
color="blue",
ls="--",
marker="o",
)
ax1[0, 1].set_title("$\\sigma_{yy}$ vs Time at $(0., 5., 0.)$")
ax1[0, 1].set_xlabel("Time / s")
ax1[0, 1].set_ylabel("$\\sigma_{yy}$ / MPa")
ax1[0, 1].legend()
ax1[0, 1].grid(True)
# Create a second y-axis for the second data series
ax2_01 = ax1[0, 1].twinx()
color = "tab:red"
ax2_01.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2_01.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2_01.tick_params(axis="y", labelcolor=color)
ax2_01.set_xlim([0, 850])
# Bottom-right plot
ax1[1, 1].plot(
stressData["Time"],
stressData["pt3_sigyy"],
#timeSteps,
#results["stress"]["pt3"].T[2] * 1e-6,
label="$\\sigma_{zz}$, TUBAF",
color="black",
ls="--",
marker="o",
)
ax1[1, 1].set_title("$\\sigma_{zz}$ vs Time $(0., 5., 0.)$")
ax1[1, 1].set_xlabel("Time / s")
ax1[1, 1].set_ylabel("$\\sigma_{zz}$ / MPa")
ax1[1, 1].legend()
ax1[1, 1].grid(True)
# Create a second y-axis for the second data series
ax2_11 = ax1[1, 1].twinx()
color = "tab:red"
ax2_11.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2_11.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2_11.tick_params(axis="y", labelcolor=color)
ax2_11.set_xlim([0, 850])
# Bottom-left plot
ax1[1, 0].plot(
stressData["Time"],
stressData["pt3_sigxy"],
#timeSteps,
#results["stress"]["pt3"].T[3] * 1e-6,
label="$\\sigma_{xy}$, TUBAF",
color="magenta",
ls="--",
marker="o",
)
ax1[1, 0].set_title("$\\sigma_{xy}$ vs Time at $(0., 5., 0.)$")
ax1[1, 0].set_xlabel("Time / s")
ax1[1, 0].set_ylabel("$\\sigma_{xy}$ / MPa")
ax1[1, 0].legend()
ax1[1, 0].grid(True)
# Create a second y-axis for the second data series
ax2_10 = ax1[1, 0].twinx()
color = "tab:red"
ax2_10.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2_10.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2_10.tick_params(axis="y", labelcolor=color)
ax2_10.set_xlim([0, 850])
# Adjust layout
plt.tight_layout()
plt.show()
Calculation of the shear and normal stresses on the fault dip¶
def stress_ns(csv_file, file, pt, theta):
stressData = pd.read_csv(csv_file)
nVector = np.array(
[math.cos(math.radians(theta)), math.sin(math.radians(theta)), 0.0]
)
# Compute the dyadic product
dyadic_product_of_the_normal_vector = np.outer(nVector, nVector)
ndyn_wp = dyadic_product_of_the_normal_vector
stress_ten = np.zeros((3, 3))
sigma_wp = np.array([0.0])
tau_wp = np.array([0.0])
sig00 = np.array([0.0])
sig11 = np.array([0.0])
sig22 = np.array([0.0])
sig01 = np.array([0.0])
sig02 = np.array([0.0])
sig12 = np.array([0.0])
stress = file[pt]
for i in range(len(stress)):
stress_ten[0, 0] = stressData[f"{pt}_sigxx"][i]
stress_ten[1, 1] = stressData[f"{pt}_sigyy"][i]
stress_ten[2, 2] = stressData[f"{pt}_sigzz"][i]
stress_ten[0, 1] = stressData[f"{pt}_sigxy"][i]
stress_ten[1, 2] = stressData[f"{pt}_sigyz"][i]
stress_ten[0, 2] = stressData[f"{pt}_sigxz"][i]
stress_ten[1, 0] = stressData[f"{pt}_sigxy"][i]
stress_ten[2, 1] = stressData[f"{pt}_sigyz"][i]
stress_ten[2, 0] = stressData[f"{pt}_sigxz"][i]
#stress_ten[0, 0] = np.array(stress)[i].T[0]
#stress_ten[1, 1] = np.array(stress)[i].T[1]
#stress_ten[2, 2] = np.array(stress)[i].T[2]
#stress_ten[0, 1] = np.array(stress)[i].T[3]
#stress_ten[0, 2] = np.array(stress)[i].T[5]
#stress_ten[1, 2] = np.array(stress)[i].T[4]
#stress_ten[1, 0] = np.array(stress)[i].T[3]
#stress_ten[2, 0] = np.array(stress)[i].T[5]
#stress_ten[2, 1] = np.array(stress)[i].T[4]
#####
s_wp = np.tensordot(stress_ten, ndyn_wp)
t_wp = np.sqrt(
np.tensordot((np.dot(stress_ten, stress_ten.transpose())), ndyn_wp)
- s_wp * s_wp
)
sigma_wp = np.append(sigma_wp, s_wp)
tau_wp = np.append(tau_wp, t_wp)
sig01 = np.append(sig01, stress_ten[0,1])
sig02 = np.append(sig02, stress_ten[0,2])
sig12 = np.append(sig12, stress_ten[1,2])
sig00 = np.append(sig00, stress_ten[0,0])
sig11 = np.append(sig11, stress_ten[1,1])
sig22 = np.append(sig22, stress_ten[2,2])
return sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp
Fault stress at P1: $(0, 0, 0)$¶
theta = 155.0
pt = "pt0"
csv_file = "stress_results_FM1.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(csv_file, file, pt, theta)
-4106470.603848967
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
timeSteps,
sig00[1:] * 1e-6,
label="$\\sigma_{xx}$, TUBAF",
color="brown",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig11[1:] * 1e-6,
label="$\\sigma_{yy}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig22[1:] * 1e-6,
label="$\\sigma_{zz}$, TUBAF",
color="black",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig01[1:] * 1e-6,
label="$\\sigma_{xy}$ TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig12[1:] * 1e-6,
label="$\\sigma_{yz}$, TUBAF",
color="blue",
ls="-",
marker="D",
)
ax1.plot(
timeSteps,
sig02[1:] * 1e-6,
label="$\\sigma_{xz}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_ylabel("$\\sigma_{\\mathrm{ij}}$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stresses at P1 $(0, 0, 0)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[17] * 1e-6,
label="Shear stress, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM1[1],
FM1[18] * 1e-6,
label="Normal stress, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau = \\sigma_n\,\\tan(\phi) + c$",
color="black",
ls="--",
marker="D",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n \quad & \quad \\tau^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P1 $(0, 0, 0)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
(0.0, 800.0)
Fault stress at P3: $(0, -0.63, -1.36)$¶
theta = 155.0
pt = "pt2"
csv_file = "stress_results_FM1.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(csv_file, file, pt, theta)
-5499905.516904995
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
timeSteps,
sig00[1:] * 1e-6,
label="$\\sigma_{xx}$, TUBAF",
color="brown",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig11[1:] * 1e-6,
label="$\\sigma_{yy}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig22[1:] * 1e-6,
label="$\\sigma_{zz}$, TUBAF",
color="black",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig01[1:] * 1e-6,
label="$\\sigma_{xy}$ TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig12[1:] * 1e-6,
label="$\\sigma_{yz}$, TUBAF",
color="blue",
ls="-",
marker="D",
)
ax1.plot(
timeSteps,
sig02[1:] * 1e-6,
label="$\\sigma_{xz}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P1 $(0, -0.63, -1.36)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[22] * 1e-6,
label="Normal Stress, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM1[1], FM1[21] * 1e-6, label="Shear Stress ENSI", color="blue", ls="--", marker="o"
) # Line plot with markers
ax1.plot(
timeSteps,
#-sig01[1:] * 1e-6,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau = \\sigma_n\,\\tan(\phi) + c$",
color="black",
ls="--",
marker="o",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P3 $(0, -0.63, -1.36)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Fault stress at P2: $(-1.5, 0., 0.)$¶
theta = 155.0
pt = "pt2"
csv_file = "stress_results_FM1.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(csv_file, file, pt, theta)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
timeSteps,
sig00[1:] * 1e-6,
label="$\\sigma_{xx}$, TUBAF",
color="brown",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig11[1:] * 1e-6,
label="$\\sigma_{yy}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig22[1:] * 1e-6,
label="$\\sigma_{zz}$, TUBAF",
color="black",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig01[1:] * 1e-6,
label="$\\sigma_{xy}$ TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
sig12[1:] * 1e-6,
label="$\\sigma_{yz}$, TUBAF",
color="blue",
ls="-",
marker="D",
)
ax1.plot(
timeSteps,
sig02[1:] * 1e-6,
label="$\\sigma_{xz}$, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P1 $(-1.5, 0., 0.)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM1[1],
FM1[20] * 1e-6,
label="Normal Stress, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM1[1], FM1[19] * 1e-6, label="Shear Stress ENSI", color="blue", ls="--", marker="o"
) # Line plot with markers
ax1.plot(
timeSteps,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau=\\sigma_n\,\\tan(\phi) + c $",
color="black",
ls="--",
marker="o",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P2 $(-1.5, 0., 0.)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Dispalcement jumps between two points (M and N) located on the left and right hand side of the fault along a line prependicular to the fault¶
def displacement_jump(file, thetap, thetam):
u_np = file["displacement"]["pt4"].T[0] * math.cos(math.radians(thetap)) + file[
"displacement"
]["pt4"].T[1] * math.sin(math.radians(thetap))
u_nm = file["displacement"]["pt5"].T[0] * math.cos(math.radians(thetam)) + file[
"displacement"
]["pt5"].T[1] * math.sin(math.radians(thetam))
u_sp = -file["displacement"]["pt4"].T[0] * math.sin(math.radians(thetap)) + file[
"displacement"
]["pt4"].T[1] * math.cos(math.radians(thetap))
u_sm = -file["displacement"]["pt5"].T[0] * math.sin(math.radians(thetam)) + file[
"displacement"
]["pt5"].T[1] * math.cos(math.radians(thetam))
return -u_np, u_nm, -u_sp, u_sm
fig, ax1 = plt.subplots(figsize=(12.5, 6.5))
thetap = 155.0
thetam = -25
u_np, u_nm, u_sp, u_sm = displacement_jump(results, thetap, thetam)
ax1.plot(
timeSteps,
u_np* 1000.0,
label="Displacement right",
ls="--",
color="blue",
marker="o",
)
ax1.plot(
timeSteps,
u_nm * 1000.0,
label="Displacement left",
ls="--",
color="red",
marker="o",
)
ax1.set_xlabel("$t$ / s")
ax1.grid(True)
ax1.set_ylabel("$u_r|_{r_\\mathrm{i}}$ / mm")
ax1.legend(loc="upper left")
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.plot(
dt.cumsum()[1:] - dt[1],
p[1:],
label="Injected pressure pressure in MPa",
color="red",
)
ax2.set_xlabel("$t-t_\\mathrm{ic}$ / s")
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.tick_params(axis="y", labelcolor=color)
Evolution of permeability $(k_{\text{d}})$ along fault dip¶
Time = [100, 453, 750]
colors = ["blue", "black", "red", "magenta"]
markers = ["d", "o"]
fig, ax = plt.subplots(figsize=(14.5, 6.5))
# Plotting loop
for k, i in enumerate(Time):
# Read permeability data
permeability_data = pvdfile_elastop_Kd_FM1.read_set_data(
i, "permeability", pointsetarray=fpoints_mid
).T[0]
# Filter out zero values
filtered_indices = permeability_data != 0
filtered_rad = rad[filtered_indices]
filtered_permeability = permeability_data[filtered_indices]
# Plot filtered data
ax.plot(
filtered_rad,
filtered_permeability,
label=f"Permeability at $t = {i}$ s (HM-EP-Kd)",
color=colors[k % len(colors)], # Cycle through colors if needed
ls="-",
marker=markers[k % len(markers)], # Cycle through markers if needed
)
# Customizations
ax.set(
xlabel="$r$ / m",
ylabel="Permeability / m$^\mathrm{2}$",
xlim=(0.0, 22.0),
)
plt.yscale("log") # Logarithmic y-axis
ax.legend()
ax.grid(True)
fig.tight_layout()
# Show the plot
plt.show()
Evolution of weak plane equivalent plastic strain along fault dip¶
Time = [420, 452, 800]
k = 0
fig, ax = plt.subplots(figsize=(14.5, 6.5))
threshold = 1e-16
for i in Time:
colors = ["blue", "black", "red", "magenta"]
markers = ["d", "o"]
peeq = pvdfile_elastop_Kd_FM1.read_set_data(
i, "EquivalentPlasticStrainWeakPlane", pointsetarray=fpoints_mid
)
peeq_data = peeq.copy()
peeq_data[peeq <= threshold] = 0.0
ax.plot(
rad[::20],
peeq_data[::20],
label="Equivalent Plastic Strain at $t = %i$ s (HM-EP-Kd)" % i,
color=colors[k],
ls="-",
marker="D",
)
k += 1
ax.set(xlim=(0.0, 22.0), autoscale_on=True)
ax.set_xlabel("$r$ / m")
ax.set_ylabel("Equivalent Plastic Strain")
ax.legend()
ax.grid(True)
fig.tight_layout()
FM2 Model¶
Import the data from PVD files¶
pvdfile_elastop_Kd_FM2 = vtuIO.PVDIO("./_out_FM2/FST_EPKd_FM2.pvd", dim=3)
Coordinates of injection, anchors, and monitoring Points¶
Anchors at poinsts at $(0, 0, 0.25)$ and $(0, 0, -0.25)$, injection at poinst at $(0, 0, 0)$, and moniroting point at $(0, -0.63, -1.36)$
anchors_point_1 = (10.0, 10.25, 0.0)
anchors_point_2 = (10.0, 9.75, 0.0)
pinjection_point = (10.0, 10.0, 0.0)
pinjection_point_left = (9.904226, 10.0446602, 0.0)
pinjection_point_right = (10.095774, 9.9553398, 0.0)
monitoring_point_p3 = (9.37, 8.64, 0.0)
monitoring_point_p3_left = (9.274225925, 8.693620825, 0.0)
monitoring_point_p3_right = (9.465774075, 8.604300455, 0.0)
monitoring_point_p2 = (10.0, 10.0, -1.5)
monitoring_point_p2_left = (9.904226, 10.0446602, -1.5)
monitoring_point_p2_right = (10.095774, 9.9553398, -1.5)
flowRate_point = {"pt0": pinjection_point}
pressure_points = {
"pt0": pinjection_point,
"pt1": monitoring_point_p3,
"pt2": monitoring_point_p2,
}
displacement_points = {
"pt0": pinjection_point,
"pt1": anchors_point_1,
"pt2": anchors_point_2,
"pt3": monitoring_point_p3,
"pt4": pinjection_point_left,
"pt5": pinjection_point_right,
"pt6": monitoring_point_p3_left,
"pt7": monitoring_point_p3_right,
"pt8": monitoring_point_p2,
"pt9": monitoring_point_p2_left,
"pt10": monitoring_point_p2_right,
}
Calling the results array of the mass flow rate, integrated pressure, displacements and stresses at the specifed points.¶
results = {}
results["MassFlowRate"] = pvdfile_elastop_Kd_FM2.read_time_series(
"MassFlowRate", pts=flowRate_point, interpolation_method="nearest"
)
results["pressure"] = pvdfile_elastop_Kd_FM2.read_time_series(
"pressure_interpolated", pts=pressure_points, interpolation_method="nearest"
)
results["displacement"] = pvdfile_elastop_Kd_FM2.read_time_series(
"displacement", pts=displacement_points, interpolation_method="nearest"
)
results["stress"] = pvdfile_elastop_Kd_FM2.read_time_series(
"sigma", pts=displacement_points, interpolation_method="nearest"
)
timeSteps = pvdfile_elastop_Kd_FM2.timesteps
Save the pressure data in a csv file¶
# Output CSV file
output_csv = "pressure_results_FM2.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 3
# Write header row: "Time" + Node IDs (e.g., pt0, pt1, pt2)
header = ["Time"] + [
f"pt{node_id}" for node_id in range(node_ids)
] # Adjust as needed
writer.writerow(header)
# Write each row: time step + scaled pressures for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Scale the pressure value by 1e-6 and append it
scaled_pressure = results["pressure"][f"pt{node_id}"][i] * 1e-6
row.append(scaled_pressure)
writer.writerow(row)
Save the mass flow rate data in a csv file¶
# Output CSV file
output_csv = "MassFlowRate_results_FM2.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 1
# Write header row: "Time" + Node IDs (e.g., pt0 ~ pt11)
header = ["Time"] + [
f"pt{node_id}" for node_id in range(node_ids)
] # Adjust as needed
writer.writerow(header)
# Write each row: time step + scaled pressures for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Scale the pressure value by 1e-6 and append it
scaled_pressure = results["MassFlowRate"][f"pt{node_id}"][i]
row.append(scaled_pressure)
writer.writerow(row)
Save the displacement data in a csv file¶
# Output CSV file
output_csv = "displacement_results_FM2.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 11 # Number of nodes (pt0 ~ pt9)
# Write header row: "Time" + Displacement Components (e.g., pt0_ux, pt0_uy, pt0_uz)
header = ["Time"]
for node_id in range(node_ids):
header.extend([f"pt{node_id}_ux", f"pt{node_id}_uy", f"pt{node_id}_uz"])
writer.writerow(header)
# Write each row: time step + displacement components for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Extract the displacement vector for the node and split into components
displacement_vector = results["displacement"][f"pt{node_id}"][i]
ux, uy, uz = (
displacement_vector[0],
displacement_vector[1],
displacement_vector[2],
)
row.extend([ux, uy, uz])
writer.writerow(row)
Save the stress data in a csv file¶
# Output CSV file
output_csv = "stress_results_FM2.csv"
# Write results to a CSV file
with open(output_csv, mode="w", newline="") as file:
writer = csv.writer(file)
node_ids = 11 # Number of nodes (pt0 ~ pt9)
# Write header row: "Time" + Displacement Components (e.g., pt0_ux, pt0_uy, pt0_uz)
header = ["Time"]
for node_id in range(node_ids):
header.extend(
[
f"pt{node_id}_sigxx",
f"pt{node_id}_sigyy",
f"pt{node_id}_sigzz",
f"pt{node_id}_sigxy",
f"pt{node_id}_sigyz",
f"pt{node_id}_sigxz",
]
)
writer.writerow(header)
# Write each row: time step + displacement components for each node
for i, time in enumerate(timeSteps):
row = [time]
for node_id in range(node_ids):
# Extract the displacement vector for the node and split into components
stress_tensor = results["stress"][f"pt{node_id}"][i]
sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_xz = (
stress_tensor[0],
stress_tensor[1],
stress_tensor[2],
stress_tensor[3],
stress_tensor[4],
stress_tensor[5],
)
row.extend([sig_xx, sig_yy, sig_zz, sig_xy, sig_yz, sig_xz])
writer.writerow(row)
Modeling results for FM1 benchmark¶
Pressure at P3: $(0, -0.63, -1.36)$¶
fig, ax = plt.subplots(figsize=(14, 6.5))
# Read the CSV file
csv_file = "pressure_results_FM2.csv"
pressureData = pd.read_csv(csv_file)
ax.plot(
FM2[1],
FM2[4] * 1e-6,
label="pressure at monitoring point, ENSI",
color="blue",
ls="--",
marker="o",
)
ax.plot(
pressureData["Time"],
pressureData["pt1"],
# instead the following two lines are used to plot the results
# timeSteps,
# results["pressure"]["pt1"] * 1e-6,
label="Pressure at monitoring Point, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Pressure at monitoring point P3 ($0.$, $-0.63$, $-1.36$)")
ax.set_xlabel("Time / s")
ax.set_ylabel("Pressure / MPa")
ax.legend()
ax.grid(True)
fig.tight_layout()
Pressure at P3: $(-1.5, 0, 0)$¶
fig, ax = plt.subplots(figsize=(14, 6.5))
ax.plot(
FM2[1],
FM2[3] * 1e-6,
label="pressure the monitoring point, ENSI",
color="blue",
ls="--",
marker="o",
)
ax.plot(
timeSteps,
results["pressure"]["pt2"] * 1e-6,
label="Pressure at the monitoring Point, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Pressure at the monitoring point ($-1.5$, $0$, $0$) s")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Pressure [MPa]")
ax.legend()
ax.grid(True)
fig.tight_layout()
Pressure along fault dip at $t= 317$¶
start_point_left = (5.32527, 0.0, 0)
end_point_left = (14.6514, 20.0, 0.0)
start_point_right = (5.34858, 0.0, 0.0)
end_point_right = (14.6748, 20.0, 0.0)
start_point_mid = (start_point_left[0] + (start_point_right[0] - start_point_left[0]) / 2., 0.0, 0.0)
end_point_mid = (end_point_left[0] + (end_point_right[0] - end_point_left[0]) / 2., 20.0, 0.0)
print(start_point_mid, end_point_mid)
stop = np.sqrt(
(end_point_mid[0] - start_point_mid[0]) ** 2.0
+ (end_point_mid[1] - start_point_mid[1]) ** 2.0
)
def divide_line(start_point, end_point, n):
fault_axis = []
for i in range(0, n + 1):
x = start_point[0] + (end_point[0] - start_point[0]) * i / n
y = start_point[1] + (end_point[1] - start_point[1]) * i / n
fault_axis.append([x, y, 0.0])
return fault_axis
npts = 100
fpoints_mid = divide_line(start_point_mid, end_point_mid, npts)
fpoints_left = divide_line(start_point_left, end_point_left, npts)
fpoints_right = divide_line(start_point_right, end_point_right, npts)
rad = np.linspace(start=0, stop=stop, num=npts + 1)
(5.336925, 0.0, 0.0) (14.6631, 20.0, 0.0)
fig, ax = plt.subplots(figsize=(14.5, 6.5))
ax.plot(
profiles_along_fault_dip_FM2[0],
profiles_along_fault_dip_FM2[1] * 1e-6,
label="pressure along the fault dip at $t=317$ s, ENSI",
color="blue",
ls="--",
marker="o",
)
itime = 317
ax.plot(
rad,
pvdfile_elastop_Kd_FM2.read_set_data(
itime, "pressure_interpolated", pointsetarray=fpoints_mid
)
* 1e-6,
label="pressure along the fault tip at $t = %i $ s, TUBAF" % itime,
color="red",
ls="-",
marker="D",
)
ax.set_title("Pressure along the fault dip at $t=317$ s")
ax.set_xlabel("Fault dip (r) / m")
ax.set_ylabel("Pressure / MPa")
ax.legend()
ax.grid(True)
fig.tight_layout()
Injection Rate¶
fig, ax = plt.subplots(figsize=(14.5, 6.5))
# Read the CSV file
csv_file = "MassFlowRate_results_FM2.csv"
MassFlowRateData = pd.read_csv(csv_file)
ax.plot(
FM2[1],
FM2[2],
label="Injection rate, ENSI",
color="blue",
ls="--",
marker="o",
)
# Initialize the total mass flow rate
total_mass_flow_rate = 0
# Loop through points
for i in range(1):
key = f"pt{i}"
total_mass_flow_rate += results["MassFlowRate"][key] * 60.0 * 2.0
MassFlowRateData["TotalMassFlowRate"] = (
MassFlowRateData.iloc[:, 1:].sum(axis=1) * 60.0 * 2.0
) # Sum across all columns except "Time"
ax.plot(
MassFlowRateData["Time"],
MassFlowRateData["TotalMassFlowRate"],
# timeSteps,
# total_mass_flow_rate,
label="Injection rate, TUBAF",
color="red",
ls="-",
marker="D",
)
ax.set_title("Injection rate versus Time")
ax.set_xlabel("Time [s]")
ax.set_ylabel("Injection Rate [L / min]")
ax.legend(loc="lower right")
ax.grid(True)
fig.tight_layout()
Differential displacements between anchors¶
fig, ax1 = plt.subplots(figsize=(14, 6.5))
csv_file = "displacement_results_FM2.csv"
displacementData = pd.read_csv(csv_file)
ax1.plot(
FM2[1],
FM2[6] * 1000.0,
label="Differential displacement ($dy$), ENSI",
color="blue",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM2[1],
FM2[7] * 1000.0,
label="Differential displacement ($dz$), ENSI",
color="red",
ls="--",
marker="o",
)
# Line plot with markers
ax1.plot(
displacementData["Time"],
(displacementData["pt1_ux"] - displacementData["pt2_ux"]) * 1e3,
# timeSteps,
# (results["displacement"]["pt0"].T[0]
# - results["displacement"]["pt1"].T[0]) * 1e3,
label="Differential displacement ($dy$), TUBAF",
color="magenta",
ls="-",
marker="D",
)
ax1.plot(
displacementData["Time"],
(displacementData["pt1_uy"] - displacementData["pt2_uy"]) * 1e3,
# timeSteps,
# (results["displacement"]["pt0"].T[1]
# - results["displacement"]["pt1"].T[1]) * 1e3,
label="Differential displacement ($dz$), TUBAF",
color="green",
ls="-",
marker="D",
)
ax1.set_title("Relative displacement of anchors")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Rotation matrix with respect to the fault dip angle¶
theta = 155
rotation_matrix = np.array(
[
[math.cos(math.radians(theta)), math.sin(math.radians(theta))],
[-math.sin(math.radians(theta)), math.cos(math.radians(theta))],
]
)
print(rotation_matrix)
[[-0.90630779 0.42261826] [-0.42261826 -0.90630779]]
Calculation of the normal and shear displacment along the fault dip¶
def displacement_ns(file, theta, pt_left, pt_right):
u_nd_left = file["displacement"][pt_left].T[0] * math.cos(
math.radians(theta)
) + file["displacement"][pt_left].T[1] * math.sin(math.radians(theta))
u_sd_left = -file["displacement"][pt_left].T[0] * math.sin(
math.radians(theta)
) + file["displacement"][pt_left].T[1] * math.cos(math.radians(theta))
##
u_nd_right = file["displacement"][pt_right].T[0] * math.cos(
math.radians(theta)
) + file["displacement"][pt_right].T[1] * math.sin(math.radians(theta))
u_sd_right = -file["displacement"][pt_right].T[0] * math.sin(
math.radians(theta)
) + file["displacement"][pt_right].T[1] * math.cos(math.radians(theta))
return abs(u_nd_left) + abs(u_nd_right), abs(u_sd_left) + abs(u_sd_right)
u_nd, u_sd = displacement_ns(results, 155, "pt1", "pt2")
fig, ax1 = plt.subplots(figsize=(14, 6.5))
ax1.plot(
FM2[1],
FM2[8] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM2[1],
FM2[10] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="-",
marker="d",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="-",
marker="d",
)
ax1.set_title("Fault displacement at $(0, 0, 0)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
u_nd, u_sd = displacement_ns(results, 155, "pt6", "pt7")
Fault displacement at P3: $(0., -0.63, -1.36)$¶
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM2[1],
FM2[14] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM2[1],
FM2[16] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_title("Fault displacement at P2 $(0., -0.63, -1.36)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
u_nd, u_sd = displacement_ns(results, 155, "pt9", "pt10")
Fault displacement at P3: $(-1.5, 0., 0.)$¶
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM2[1],
FM2[11] * 1000.0,
label="Shear displacement, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM2[1],
FM2[13] * 1000.0,
label="Normal displacement, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
u_sd * 1000.0,
label="Shear displacement, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
u_nd * 1000.0,
label="Normal displacement, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.set_title("Fault displacement at P2 $(-1.5, 0., 0.)$")
ax1.set_xlabel("Time")
ax1.set_ylabel("Displacement [mm]")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Calculation of the shear and normal stresses on the fault dip¶
def stress_ns(csv_file, file, pt, theta):
stressData = pd.read_csv(csv_file)
nVector = np.array(
[math.cos(math.radians(theta)), math.sin(math.radians(theta)), 0.0]
)
# Compute the dyadic product
dyadic_product_of_the_normal_vector = np.outer(nVector, nVector)
ndyn_wp = dyadic_product_of_the_normal_vector
stress_ten = np.zeros((3, 3))
sigma_wp = np.array([0.0])
tau_wp = np.array([0.0])
sig00 = np.array([0.0])
sig11 = np.array([0.0])
sig22 = np.array([0.0])
sig01 = np.array([0.0])
sig02 = np.array([0.0])
sig12 = np.array([0.0])
stress = file[pt]
for i in range(len(stress)):
stress_ten[0, 0] = stressData[f"{pt}_sigxx"][i]
stress_ten[1, 1] = stressData[f"{pt}_sigyy"][i]
stress_ten[2, 2] = stressData[f"{pt}_sigzz"][i]
stress_ten[0, 1] = stressData[f"{pt}_sigxy"][i]
stress_ten[1, 2] = stressData[f"{pt}_sigyz"][i]
stress_ten[0, 2] = stressData[f"{pt}_sigxz"][i]
stress_ten[1, 0] = stressData[f"{pt}_sigxy"][i]
stress_ten[2, 1] = stressData[f"{pt}_sigyz"][i]
stress_ten[2, 0] = stressData[f"{pt}_sigxz"][i]
# stress_ten[0, 0] = np.array(stress)[i].T[0]
# stress_ten[1, 1] = np.array(stress)[i].T[1]
# stress_ten[2, 2] = np.array(stress)[i].T[2]
# stress_ten[0, 1] = np.array(stress)[i].T[3]
# stress_ten[0, 2] = np.array(stress)[i].T[5]
# stress_ten[1, 2] = np.array(stress)[i].T[4]
# stress_ten[1, 0] = np.array(stress)[i].T[3]
# stress_ten[2, 0] = np.array(stress)[i].T[5]
# stress_ten[2, 1] = np.array(stress)[i].T[4]
#####
s_wp = np.tensordot(stress_ten, ndyn_wp)
t_wp = np.sqrt(
np.tensordot((np.dot(stress_ten, stress_ten.transpose())), ndyn_wp)
- s_wp * s_wp
)
sigma_wp = np.append(sigma_wp, s_wp)
tau_wp = np.append(tau_wp, t_wp)
sig01 = np.append(sig01, stress_ten[0, 1])
sig02 = np.append(sig02, stress_ten[0, 2])
sig12 = np.append(sig12, stress_ten[1, 2])
sig00 = np.append(sig00, stress_ten[0, 0])
sig11 = np.append(sig11, stress_ten[1, 1])
sig22 = np.append(sig22, stress_ten[2, 2])
return sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp
Fault stress at P1: $(0, 0, 0)$¶
theta = 155.0
pt = "pt0"
csv_file = "stress_results_FM2.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(
csv_file, file, pt, theta
)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM2[1],
FM2[17] * 1e-6,
label="Shear stress, ENSI",
color="blue",
ls="--",
marker="o",
)
ax1.plot(
FM2[1],
FM2[18] * 1e-6,
label="Normal stress, ENSI",
color="red",
ls="--",
marker="o",
)
ax1.plot(
timeSteps,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau = \\sigma_n\,\\tan(\phi) + c$",
color="black",
ls="--",
marker="D",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n \quad & \quad \\tau^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P1 $(0, 0, 0)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
(0.0, 800.0)
Fault stress at P3: $(0, -0.63, -1.36)$¶
theta = 155.0
pt = "pt3"
csv_file = "stress_results_FM2.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(
csv_file, file, pt, theta
)
fig, ax1 = plt.subplots(figsize=(14, 6.5))
ax1.plot(
FM2[1],
FM2[22] * 1e-6,
label="Normal Stress, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM2[1], FM2[21] * 1e-6, label="Shear Stress ENSI", color="blue", ls="--", marker="o"
) # Line plot with markers
ax1.plot(
timeSteps,
# -sig01[1:] * 1e-6,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau = \\sigma_n\,\\tan(\phi) + c$",
color="black",
ls="--",
marker="o",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P3 $(0, -0.63, -1.36)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Fault stress at P2: $(-1.5, 0., 0.)$¶
theta = 155.0
pt = "pt8"
csv_file = "stress_results_FM2.csv"
file = results["stress"]
sig01, sig02, sig12, sig00, sig11, sig22, sigma_wp, tau_wp = stress_ns(
csv_file, file, pt, theta
)
fig, ax1 = plt.subplots(figsize=(13.5, 7))
ax1.plot(
FM2[1],
FM2[20] * 1e-6,
label="Normal Stress, ENSI",
color="red",
ls="--",
marker="o",
) # Line plot with markers
ax1.plot(
FM2[1], FM2[19] * 1e-6, label="Shear Stress ENSI", color="blue", ls="--", marker="o"
) # Line plot with markers
ax1.plot(
timeSteps,
tau_wp[1:] * 1e-6,
label="Shear stress, TUBAF",
color="magenta",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6,
label="Normal stress, TUBAF",
color="green",
ls="--",
marker="D",
)
ax1.plot(
timeSteps,
-sigma_wp[1:] * 1e-6 * math.tan(math.radians(22)),
label="$\\tau=\\sigma_n\,\\tan(\phi) + c $",
color="black",
ls="--",
marker="o",
)
ax1.set_ylabel("$\\sigma^{\\mathrm{wp}}_n$ / MPa")
ax1.set_xlabel("Time / s")
ax1.set_title("Fault stress at P2 $(-1.5, 0., 0.)$")
ax1.legend()
ax1.grid(True)
# Create a second y-axis for the second data series
ax2 = ax1.twinx()
color = "tab:red"
ax2.set_ylabel("Injection Pressure [MPa]", color=color) # For the second data series
ax2.plot(injectionPressure[0], injectionPressure[4] / 1e6, color=color)
ax2.tick_params(axis="y", labelcolor=color)
ax2.set_xlim([0, 800])
# Set x-axis limits
(0.0, 800.0)
Evolution of permeability $(k_{\text{d}})$ along fault dip¶
Time = [100, 453, 800]
k = 0
fig, ax = plt.subplots(figsize=(14.5, 6.5))
for i in Time:
colors = ["blue", "black", "red", "magenta"]
markers = ["d", "o"],
ax.plot(
rad,
pvdfile_elastop_Kd_FM2.read_set_data(
i, "permeability", pointsetarray=fpoints_mid
).T[0],
label="Permeability at $t = %i$ s (HM-EP-Kd)" % i,
color=colors[k],
ls="-",
marker="D",
)
k += 1
ax.set(xlim=(0.0, 22.0), autoscale_on=True)
plt.yscale("log") # Set the y-axis to logarithmic scale
ax.set_xlabel("$r$ / m")
ax.set_ylabel("Permeability / m$^\mathrm{2}$")
ax.legend()
ax.grid(True)
fig.tight_layout()
Time = [100, 453, 800]
colors = ["blue", "black", "red", "magenta"]
markers = ["d", "o"]
fig, ax = plt.subplots(figsize=(14.5, 6.5))
# Plotting loop
for k, i in enumerate(Time):
# Read permeability data
permeability_data = pvdfile_elastop_Kd_FM2.read_set_data(
i, "permeability", pointsetarray=fpoints_mid
).T[0]
# Filter out zero values
filtered_indices = permeability_data != 0
filtered_rad = rad[filtered_indices]
filtered_permeability = permeability_data[filtered_indices]
# Plot filtered data
ax.plot(
filtered_rad,
filtered_permeability,
label=f"Permeability at $t = {i}$ s (HM-EP-Kd)",
color=colors[k % len(colors)], # Cycle through colors if needed
ls="-",
marker=markers[k % len(markers)], # Cycle through markers if needed
)
# Customizations
ax.set(
xlabel="$r$ / m",
ylabel="Permeability / m$^\mathrm{2}$",
xlim=(0.0, 22.0),
)
plt.yscale("log") # Logarithmic y-axis
ax.legend()
ax.grid(True)
fig.tight_layout()
# Show the plot
plt.show()
Discussion of results for both fault model benchmarks¶
In this study, the results are analyzed in terms of pressure variations at the monitoring points P2 and P3 and compared with ENSI's findings. According to the definition of FM1, the pressure at the monitoring points remains constant in the absence of shear failure and begins to evolve only upon the onset of dilation within the fault. In contrast, under the FM2 framework, pressure evolution occurs from the initiation of well pressure injection. In FM1, an induced aperture change of $28 \mu \text{m}$ enhances the permeability of the fault as soon as shear or tensile failure is initiated, leading to abrupt pressure variations at the peak injection pressure of $6.3$ MPa. Furthermore, this induced aperture change results in increased permeability in the FM1 benchmark, causing the maximum pressure at $t=420.5$ s to exceed that observed in the FM2 benchmark. Moreover, the comparison of pressure evolution along the fault direction at $t=317$ s demonstrates a reasonable agreement with ENSI's findings, and conforms to the concepts of the fault models.
When comparing the injection flow rates, in the case of FM2, the injection flow rate increases with each stepwise increase in injection pressure at the injection points. This trend results from the increasing pressure gradient and permeability induced by fracture opening. The maximum injection flow rate is reached at the peak injection pressure of $6.3$ MPa, coinciding with the onset of dilation. Subsequently, an abrupt decrease occurs at $t=453.5$ s when the injection pressure drops to $3.382$ MPa.
In contrast, for FM1, since neither normal nor shear opening occurs before the onset of dilation, the injection flow rate remains unchanged compared to the FM2 benchmark. However, the maximum injection flow rate is attained at the peak injection pressure of $6.3$ MPa at $t=420.5$ s, indicating the initiation of dilation, accompanied by a sudden change in the permeability of the fault. A subsequent decrease follows at $t=453.5$ s as the injection pressure drops to $3.382$ MPa, ultimately leading to the hydraulic closure of the discontinuity.
In terms of the mechanical results, the relative displacements at the anchors are plotted. Furthermore, the normal and shear displacements, as well as the corresponding stresses, are also presented. The comparison of the results with ENSI's findings demonstrates an acceptable agreement between the normal displacements and stresses at the faults.
Moreover, the evolution of permeability along the fault direction is presented. According to the FM1 concept, an abrupt change in permeability occurs when shearing initiates at approximately $t=420.5$ s, coinciding with the peak injection pressure. In contrast, in the FM2 benchmark, permeability evolves from the onset of pressure injection, with normal opening, which is then accompanied by shear opening, leading to permeability enhancement.
Conclusion¶
This work focuses on the development and adaptation of numerical simulations based on two different fault models. Accordingly, two distinct parameter sets are applied to account for shear-induced hydraulic aperture changes and the resulting variations in fracture transmissivity.
FM1: In this model, the fault undergoes elastic (normal) and plastic (shear) opening only when shear failure occurs. Initially, the normal deformation remains within the elastic regime, meaning that any changes in normal stress result in reversible displacements without permanent deformation. However, shear-induced opening follows a plastic response, which is activated only when the applied shear stress exceeds the material’s failure threshold. This behavior reflects the fault's mechanical response to stress redistribution, where normal loading alone does not contribute to permanent dilation. Instead, significant shear displacement and associated fault opening are only triggered once the fault reaches a critical stress state. This approach effectively represents fault activation mechanics in fractured rock masses, particularly under shear-driven deformation scenarios in geomechanical and hydro-mechanical coupled systems. This model is developed based on insights from modeling fault activation experiments previously carried out at the Tournemire site in Southern France.
FM2: In this mode, the fracture initially possesses a pre-existing aperture, allowing for elastic (normal) opening as pore pressure propagates. In contrast, shear-induced (dilatant) opening occurs only upon the onset of shear failure. This distinction highlights the fundamental difference between normal and shear deformation mechanisms in fractured rock masses. While normal opening is primarily governed by fluid pressure variations, leading to reversible deformation, shear dilation is controlled by the mobilization of shear stress surpassing the failure threshold, resulting in irreversible structural changes. The interaction between mechanical deformation and pore pressure evolution significantly influences fracture permeability, fault reactivation, and hydro-mechanical coupling in subsurface environments. Understanding these processes is essential for evaluating fluid flow in fractured reservoirs, assessing fault stability in geotechnical applications, and predicting induced seismicity in subsurface engineering projects.
In this study, an embedded weakness plane model is developed based on the Coulomb failure criterion, implemented in MFront, and extended to the hydro-mechanical (HM) process in OpenGeoSys. This includes an embedded fracture permeability model formulated based on the cubic law for fracture flow, which accounts for normal and shear opening as well as the induced hydraulic aperture. The results are presented in terms of pressure variation over time at monitoring points, injection flow rate, relative displacements at anchors, normal and shear displacements and stresses at monitoring points, and the evolution of permeability along the fault direction for both fault model benchmarks."
In summary, the comparison of the current results with ENSI's findings confirms the applicability of the extended elastoplastic weakness plane model, implemented in MFront, within the hydromechanical framework of the multi-physics solver OpenGeoSys. The model can be considered suitable for investigating hydro-mechanical fault reactivation at a repository site in argillaceous claystone. However, experimental investigation and subsequent validation are necessary to ensure its accuracy.
References¶
[1] Rutqvist, J., Graupner, B., Guglielmi, Y., Kim, T., Maßmann, J., Nguyen, T.S., Park, J.W., Shiu, W., Urpi, L., Yoon, J.S. and Ziefle, G., 2020. An international model comparison study of controlled fault activation experiments in argillaceous claystone at the Mont Terri Laboratory. International Journal of Rock Mechanics and Mining Sciences, 136, p.104505.
[2] Seyedi, D.M., Plúa, C., Vitel, M., Armand, G., Rutqvist, J., Birkholzer, J., Xu, H., Guo, R., Thatcher, K.E., Bond, A.E. and Wang, W., 2021. Upscaling THM modeling from small-scale to full-scale in-situ experiments in the Callovo-Oxfordian claystone. International Journal of Rock Mechanics and Mining Sciences, 144, p.104582.
[3] Ghabezloo, S. and Sulem, J., 2009. Stress dependent thermal pressurization of a fluid-saturated rock. Rock Mechanics and Rock Engineering, 42, pp.1-24.
[4] Rutqvist, J., Zheng, L., Chen, F., Liu, H. H., & Birkholzer, J. 2014. Modeling of coupled thermo-hydro-mechanical processes with links to geochemistry associated with bentonite-backfilled repository tunnels in clay formations. Rock Mechanics and Rock Engineering, 47, 167-186.
[5] Davies, R., Foulger, G., Bindley, A. and Styles, P., 2013. Induced seismicity and hydraulic fracturing for the recovery of hydrocarbons. Marine and petroleum geology, 45, pp.171-185.
[6] Flewelling, S.A., Tymchak, M.P. and Warpinski, N., 2013. Hydraulic fracture height limits and fault interactions in tight oil and gas formations. Geophysical Research Letters, 40(14), pp.3602-3606.
[7] Lei, X., Huang, D., Su, J., Jiang, G., Wang, X., Wang, H., Guo, X. and Fu, H., 2017. Fault reactivation and earthquakes with magnitudes of up to Mw4.7 induced by shale-gas hydraulic fracturing in Sichuan Basin, China. Scientific reports, 7(1), p.7971.
[8] Zill, F., Lüdeling, C., Kolditz, O., & Nagel, T. 2021. Hydro-mechanical continuum modelling of fluid percolation through rock salt. International Journal of Rock Mechanics and Mining Sciences, 147, 104879.
[9] Urpi, L., Graupner, B., Wang, W., Nagel, T., & Rinaldi, A. P. 2020. Hydro-mechanical fault reactivation modeling based on elasto-plasticity with embedded weakness planes. Journal of Rock Mechanics and Geotechnical Engineering, 12(4), 877-885.
[10] Corkum, A.G. and Martin, C.D., 2007. Modelling a mine-by test at the Mont Terri rock laboratory, Switzerland. International Journal of Rock Mechanics and Mining Sciences, 44(6), pp.846-859.
[11] Guglielmi, Y., Elsworth, D., Cappa, F., Henry, P., Gout, C., Dick, P. and Durand, J., 2015. In situ observations on the coupling between hydraulic diffusivity and displacements during fault reactivation in shales. Journal of Geophysical Research: Solid Earth, 120(11), pp.7729-7748.
[12] Singh, B. (1973, July). Continuum characterization of jointed rock masses: Part I—The constitutive equations. In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts (Vol. 10, No. 4, pp. 311-335). Pergamon.
[13] Chang, L. and Konietzky, H., 2018. Application of the Mohr-Coulomb yield criterion for rocks with multiple joint sets using fast Lagrangian analysis of continua 2D (FLAC2D) software. Energies, 11(3), p.614.
[14] Davis, R.O. and Selvadurai, A.P., 2005. Plasticity and geomechanics. Cambridge university press.
[15] Bilke, L., Naumov, D., Wang, W., Fischer, T., Kiszkurno, F.K., Lehmann, C., Jäschke, M., Zill, F., Buchwald, J., Grunwald, N., Kessler, K., Aubry, L., Dörnbrack, M., Nagel, T., Ahrendt, L., Kaiser, S., & Meisel, T. (2025) OpenGeoSys [Software], version 6.5.4, Zenodo.
[16] Helfer, T., Michel, B., Proix, J.M., Salvo, M., Sercombe, J. and Casella, M., 2015. Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. Computers & Mathematics with Applications, 70(5), pp.994-1023.
[17] Helfer, T., Bleyer, J., Frondelius, T., Yashchuk, I., Nagel, T. and Naumov, D., 2020. The MFrontGenericInterfaceSupport project. Journal of Open Source Software, 5(48), pp.1-8.